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Abstract 

The phase space flow of a dynamical system, leading to the solution of Linear Programming (LP) 
problems, is explored as an example of complexity analysis in an analog computation framework. In 
this framework, computation by physical devices and natural systems, evolving in continuous phase 
space and time (in contrast to the digital computer where these are discrete) , is explored. A Gaussian 
ensemble of LP problems is studied. The convergence time of a flow to the fixed point representing 
the optimal solution, is computed. The cumulative distribution function of the convergence time is 
calculated in the framework of Random Matrix Theory (RMT) in the asymptotic limit of large problem 
size. It is found to be a scaling function, of the form obtained in the theories of critical phenomena 
and Anderson localization. It demonstrates a correspondence between problems of Computer Science 
and Physics. 

PACS numbers: 5.45-a, 89.79+c, 89.75.D 

An analog computer is a physical device that performs computation, evolving in con- 
tinuous time and phase space; its evolution in phase space can be modeled by dynamical 
systems (DS) P, the way classical systems such as particles moving in a potential (or 
electric circuits), are modeled. This description makes a large set of analytical tools and 
physical intuition, developed for dynamical systems, applicable to the analysis of analog 
computers. In contrast, the evolution of a digital computer is described by a dynamical 
system, discrete both in its phase space and in time. The most relevant examples of analog 
computers are VLSI devices implementing neural networks j2j, or neuromorphic systems 
0, whose structure is directly motivated by the workings of the brain. Various processes 
taking place in living cells can be considered as analog computation |4j. Dynamical systems 
(described by ordinary differential equations) are also used to solve computational prob- 
lems El El ■ The standard theory of computation and computational complexity |H1 deals 
with computation in discrete time and phase space, and is inadequate for the description 
of such systems. For the analysis of computation by analog devices a theory that is valid 
in continuous time and phase space is required. Since the systems in question are physical 



1 



systems, the computation time is the time required for a system to reach the vicinity of 
an attractor (a stable fixed point in the present work) combined with the time required 
to verify that it indeed reached this vicinity. This time is the elapsed time measured by a 
clock, contrary to standard computation theory, where it is the number of steps. 

In the exploration of physical systems, it is sometimes much easier to study statistical 
ensembles of systems, estimating their typical behavior using statistical methods [Ul llUlITT] . 
Ensembles of systems modeling the dynamics of populations were studied as well jl2( ll.Sj. 
The statistical theories describe many general features of the problems that are inversti- 
gated, but specific systems require special attention j^EJ. In this letter a statistical 
theory is used to calculate the computational complexity of a standard representative prob- 
lem, namely Linear Programming (LP), as solved by a DS. A detailed version was published 
in the computer science literature fl4j . 

In two recent papers we have proposed a framework for computing with DS that converge 
exponentially to fixed points ^Hl- For such systems it is natural to consider the attracting 
fixed point as the output. The input can be modeled in various ways. One possible choice 
is the initial condition. This is appropriate when the aim of the computation is to decide 
to which attractor out of many possible ones the system flows ^Hl- Here, as in ^Sj, the 
parameters on which the system of DS depends (e.g., the parameters appearing in the vector 
field F in (^) are the input. 

The basic entity of the computational model is a dynamical system that may be 
defined by a set of Ordinary Differential Equations (ODEs) 



where x is an n-dimensional vector, and F is an n-dimensional smooth vector field, which 
converges exponentially to a fixed point. Eq. solves a computational problem as follows: 
Given an instance of the problem, the parameters of the vector field F are set, and it is 
started from some pre-determined initial condition. The result of the computation is then 
deduced from the fixed point that the system approaches. 

In our model we assume we have a physical implementation of the flow equation (^. 
Thus, the vector field F need not be computed, and the computation time is determined by 
the convergence time to the attractive fixed point. In other words, the time of flow to the 
vicinity of the attractor is a good measure of complexity, namely the computational effort, 
for the class of continuous dynamical systems introduced above [T3] . 

In this letter we will consider real inputs, as the ones found in physical experiments, 
and that are studied in the BSS model ^7j. For computational models defined on the real 
numbers, worst case behavior, that is traditionally studied in computer science, can be ill 
defined and lead to infinite computation times, in particular, for some methods for solving 
LP ^Zl^j. Therefore, we compute the distribution of computation times for a probabilistic 
model of LP instances with Gaussian distribution of the data like in |191 l2Uj . Ill-defined 
instances constitute a set of zero measure in our probability ensemble, and need not be 
concerned about. 

The computational complexity of the method presented here is 0{n\ogn), compared to 
©(n^'^logn) found for standard interior point methods [211. The basic reason is that for 
standard methods (such as interior point methods), the major component of the complexity 
of each iteration is 0{n^) due to matrix decomposition and inversion of the constraint 
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matrix, while here, because of its analog nature, the system just flows according to its 
equations of motion (which need not be computed). 

Since we consider the evolution of a vector filed, our model is inherently parallel. There- 
fore, to make the analog vs. digital comparison entirely fair, we should compare the com- 
plexity of our method to that of the best parallel algorithm. The latter can reduce the 
0{n^) time needed for matrix decomposition/inversion to poly logarithmic time (for well- 
posed problems), at the cost of 0{n'^'^) processors [22) while our system of equations ^ 
uses only 0{n) variables. 

Linear programming is a P-complete problem [H], i.e. it is representative of all problems 
that can be solved in polynomial time. The standard form of LP is to find 

max{c^x : X eWC',Ax = b,x>0} (2) 

where c € M", b G R*", A £ R™^" and m < n. The set generated by the constraints in Q 
is a polyheder. If a bounded optimal solution exists, it is obtained at one of its vertices. 
The vector defining this optimal vertex can be decomposed (in an appropriate basis) in the 
form X = {xj\f, xq) where xj\f = is an n — m component vector, while xq = B~^b > is 
an m component vector, and B is the m x m matrix whose columns are the columns of A 
with indices identical to the ones of xs- Similarly, we decompose A = (N,B). 

A flow of the form (pQ) converging to the optimal vertex, introduced by Faybusovich [Hj 
will be studied here. Its vector fleld F is a projection of the gradient of the cost function 
c^x onto the constraint set, relative to a Riemannian metric which enforces the positivity 
constraints x > 6,. It is given by 

F{x) = [X - XA^{AXA^)-^AX] c , (3) 

where X is the diagonal matrix Diag(xi . . . x„). The nm + n entries of A and c, namely, the 
parameters of the vector fleld constitute the input; as in other models of computation, we 
ignore the time it takes to "load" the input, since this step does not reflect the complexity of 
the computation being performed, either in analog or digital computation. It was shown in 
[23j that the flow equations given by (pQ) and (jSJ are, in fact, part of a system of Hamiltonian 
equations of motion of a completely integrable system of a Toda type. Therefore, like the 
Toda system, it is integrable with the formal solution |Bj 

x,{t) = xM exp I -A,t + a,, log ^l±l^zllim (4) 

\ j^l Xj-^-n-m{^) I 

[i = 1, . . . ,n — m), that describes the time evolution of the n — m independent variables 
xj\f{t), in terms of the variables 2;g(t). In Q Xj(0) and Xj+n-miO) are components of the 
initial condition, Xj+n-m{t) are the xb components of the solution, aji = —{B^^N)ji is an 
m X (n — m) matrix, while 

m 

Ai = -Ci - ^ Cjaji ■ (5) 
i=i 

For the decomposition x = {xj\f, xg) used for the optimal vertex Aj>0 i = 1, . . . ,n — m , 
and xj\f{t) converges to 0, while XQ^t) converges to x* = B~^b. Note that the analytical 
solution is only a formal one, and does not provide an answer to the LP instance, since 
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the Aj depend on the partition of A, and only relative to a partition corresponding to a 
maximum vertex are all the Aj positive. 

The second term in Q), when it is positive, is a kind of "barrier": Aji must be larger 
than the barrier before Xi can decrease to zero. In the following we ignore the contribution 
of the initial condition and denote the value of this term in the infinite time limit by 

m 

A = Zl«iilog^i+n-m- (6) 

Note that although one of the may vanish, in the probabilistic ensemble studied here 
such an event is of measure zero and therefore should not be considered. In order for x{t) 
to be close to the maximum vertex we must have Xi{t) < e for i = 1, . . . ,n — m for some 
small positive e, namely exp(— Ajt + < e , for i = 1, . . . ,n — m. Therefore we consider 

/ A |loge|\ 

T = max — + — — , (7) 



Ai A^ 
as the computation time. We denote 

Amin = minAj, /3max = max/3j . (8) 

The Aj can be arbitrarily small when the inputs are real numbers, but in the probabilistic 
model, "bad" instances are rare as is clear from H1U|) . 

The ensemble we analyze consists of LP problems in which the components of (A, 6, c) 
are independent identically distributed (i.i.d.) random variables taken from the standard 
Gaussian distribution with mean and unit variance. With the introduction of a prob- 
abilistic model of LP instances, Amin, /3max and T become random variables. Since the 
expression for Aj, equation ©, is independent of 6, its distribution is independent of h. For 
a given realization of A and c, with a partition of A into (A^, 5) such that Aj > 0, there 
exists a vector h such that the resulting polyheder has a bounded optimal solution. Since 
h in our probabilistic model is independent of A we obtain: 

'P(Ajnm < A|Amin > 0, LP instance has a bounded maximum vertex) = 'P(Amin < A|Amin > 
0). 

We wish to compute the probability distribution of Amin for instances with a bounded 
solution, when Amin > 0, denoted by Vy^min > ^\^min > 0). It turns out that it is much 
easier to analytically calculate the probability distribution of Amin for a given partition of 
the matrix A. In the probabilistic model we defined, 'P(Ajnin > A| Amin > 0) is proportional 
to the probability that Amin > A for a fixed partition @. Let the index 1 stand for the 
partition where B is taken from the last m columns of A. In Jl] we proved, using the 
symmetry resulting from the identity of the Gaussian variables, that 

P(Am.n > AlAmin > 0) = ^(^^ ^ > ^) (Q) 

for A > 0, where Ammi = min{Aj | Aj are computed relative to the partition 1} and 
P(A™„i > 0) = 1/2— . 

Integrating over the Gaussian variables of the ensemble, the probability P(Ammi > A) 
was computed in for a specific partition of A in the large (n, m) asymptotic limit, 
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Figure 1: T{x^) is plotted against the variable xa, for n = 2m. There is very good agree- 
ment with the analytical results, improving as m increases, as expected for an asymptotic 
result. 



making use of methods of random matrix theory. Given V{/S.mini > then P(Ajnin > 
A|Ajnin > 0) is obtained with the help of Q. In the large (n,m) limit the probability 
P(A„i„ < A|A™„ > 0) = :f("''")(A) is of the scaling form 

^K'")(A) = l-e^Aerfc(xA) = H^a)- (10) 

with the scaling variable x/\{n,m) = — l) y/mlS.. The scaling function T contains 

all asymptotic information on A. The distribution J^{xa) is very wide and does not have 
a mean. Also the average of diverges. 

In order to the demonstrate this result numerically, we generated LP instances where 
{A, b, c) are random Gaussian variables and solved for them the LP problem with the IMSL 
C hbrary. We obtained an estimate of .^("'''"^(A) = V{Amin < A|Amm > 0), and of the 
corresponding cumulative distribution functions of the barrier /3max and T. In Fig. 1 the 
numerical results are compared with the analytical formula H10|) . 

The existence of scaling functions like for the barrier Pmax, that is the maximum 
of the Pi defined by © and for T defined by was verified numerically (see Fig. 2 for 
T). In particular for fixed m/n, we found that 

v{i/p^a. < = :f[%ZI (1//3) = :F,/f,_^ (x^) (11) 

and 

V{l/T < 1/t) ^ Ti'jfil/t) = TyrixT). (12) 

The scaling variables are Xj3 ~ m/(3 and xt ~ mlogm/t. The scaling functions 
and p2() imply the asymptotic behavior 

1/Ainin ~ /3max ~ "I, T ~ m log m (13) 

with "high probability" [H]. 
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Figure 2: Tiix{xt) as a function of the variable xt = mlogm/t for the same instances as 
Fig. 1. 

In this letter we computed the problem size dependence of the distributions of quantities 
that govern the convergence of a DS that solves the LP problem j^j. To the best of our 
knowledge, this is the first time such distributions are computed. In particular, knowledge 
of the distribution functions enables to obtain the "high probability" behavior ()13p. and 
the moments (if these exist). The main result of the present work is that the distribution 
functions of the convergence rate, Amm; the barrier fimax and the computation time T are 
scaling functions; i.e., in the asymptotic limit of large (n, m), each depends on the problem 
size only through a scaling variable. In other words these are not arbitrary functions of the 
three variables, but each is a function only of one variable, xa, a^/j or xt- The distribution 
function of A^m was calculated analytically, and the result was verified numerically. The 
scaling functions, even if known only numerically, can be useful for the understanding of 
the behavior for large values of (n, m) that are beyond the limits of numerical simulations. 

In this letter the distribution functions of various quantities that characterize the com- 
putational complexity, were found to be scaling functions in the large (n, m) limit. This is 
analogous to the situation found for the central limit theorem, for critical phenomena [21] 
and for Anderson localization 25,, in spite of the very different nature of these problems. 
It is demonstrated here how for the implementation of the LP problem on a physical de- 
vice, methods used in theoretical physics enable to describe the distribution of computation 
times in a simple and physically transparent form. Based on our experience with certain 
universality properties of rectangular and chiral random matrix models 2^ , we expect some 
universality for computational problems, that should be explored. The obvious questions 
are: Is the Gaussian nature of the ensemble unimportant in analogy with [26|? Are there 
universality classes |24| of analog computational problems, and if they exist, what are they? 
Are these analogous to the classification of jl2|? We believe it can be instructive to explore 
computational problems using methodologies of theoretical physics as was demonstrated 
here for linear programming. 
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